<!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html>
 <head>
  <title>sggglm.f</title>
 <meta name="generator" content="emacs 21.3.1; htmlfontify 0.20">
<style type="text/css"><!-- 
body { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default   { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.default a { background: rgb(255, 255, 255);  color: rgb(0, 0, 0);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.string   { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.string a { color: rgb(188, 143, 143);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
span.comment   { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: none; }
span.comment a { color: rgb(178, 34, 34);  background: rgb(255, 255, 255);  font-style: normal;  font-weight: 500;  font-stretch: normal;  font-family: adobe-courier;  font-size: 11pt;  text-decoration: underline; }
 --></style>

 </head>
  <body>

<pre>
      SUBROUTINE <a name="SGGGLM.1"></a><a href="sggglm.f.html#SGGGLM.1">SGGGLM</a>( N, M, P, A, LDA, B, LDB, D, X, Y, WORK, LWORK,
     $                   INFO )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  -- LAPACK driver routine (version 3.1) --
</span><span class="comment">*</span><span class="comment">     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
</span><span class="comment">*</span><span class="comment">     November 2006
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Scalar Arguments ..
</span>      INTEGER            INFO, LDA, LDB, LWORK, M, N, P
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Array Arguments ..
</span>      REAL               A( LDA, * ), B( LDB, * ), D( * ), WORK( * ),
     $                   X( * ), Y( * )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Purpose
</span><span class="comment">*</span><span class="comment">  =======
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  <a name="SGGGLM.19"></a><a href="sggglm.f.html#SGGGLM.1">SGGGLM</a> solves a general Gauss-Markov linear model (GLM) problem:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          minimize || y ||_2   subject to   d = A*x + B*y
</span><span class="comment">*</span><span class="comment">              x
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where A is an N-by-M matrix, B is an N-by-P matrix, and d is a
</span><span class="comment">*</span><span class="comment">  given N-vector. It is assumed that M &lt;= N &lt;= M+P, and
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">             rank(A) = M    and    rank( A B ) = N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Under these assumptions, the constrained equation is always
</span><span class="comment">*</span><span class="comment">  consistent, and there is a unique solution x and a minimal 2-norm
</span><span class="comment">*</span><span class="comment">  solution y, which is obtained using a generalized QR factorization
</span><span class="comment">*</span><span class="comment">  of the matrices (A, B) given by
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     A = Q*(R),   B = Q*T*Z.
</span><span class="comment">*</span><span class="comment">           (0)
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  In particular, if matrix B is square nonsingular, then the problem
</span><span class="comment">*</span><span class="comment">  GLM is equivalent to the following weighted linear least squares
</span><span class="comment">*</span><span class="comment">  problem
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">               minimize || inv(B)*(d-A*x) ||_2
</span><span class="comment">*</span><span class="comment">                   x
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  where inv(B) denotes the inverse of B.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  Arguments
</span><span class="comment">*</span><span class="comment">  =========
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  N       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of rows of the matrices A and B.  N &gt;= 0.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  M       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of columns of the matrix A.  0 &lt;= M &lt;= N.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  P       (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The number of columns of the matrix B.  P &gt;= N-M.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  A       (input/output) REAL array, dimension (LDA,M)
</span><span class="comment">*</span><span class="comment">          On entry, the N-by-M matrix A.
</span><span class="comment">*</span><span class="comment">          On exit, the upper triangular part of the array A contains
</span><span class="comment">*</span><span class="comment">          the M-by-M upper triangular matrix R.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDA     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array A. LDA &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  B       (input/output) REAL array, dimension (LDB,P)
</span><span class="comment">*</span><span class="comment">          On entry, the N-by-P matrix B.
</span><span class="comment">*</span><span class="comment">          On exit, if N &lt;= P, the upper triangle of the subarray
</span><span class="comment">*</span><span class="comment">          B(1:N,P-N+1:P) contains the N-by-N upper triangular matrix T;
</span><span class="comment">*</span><span class="comment">          if N &gt; P, the elements on and above the (N-P)th subdiagonal
</span><span class="comment">*</span><span class="comment">          contain the N-by-P upper trapezoidal matrix T.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LDB     (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The leading dimension of the array B. LDB &gt;= max(1,N).
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  D       (input/output) REAL array, dimension (N)
</span><span class="comment">*</span><span class="comment">          On entry, D is the left hand side of the GLM equation.
</span><span class="comment">*</span><span class="comment">          On exit, D is destroyed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  X       (output) REAL array, dimension (M)
</span><span class="comment">*</span><span class="comment">  Y       (output) REAL array, dimension (P)
</span><span class="comment">*</span><span class="comment">          On exit, X and Y are the solutions of the GLM problem.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  WORK    (workspace/output) REAL array, dimension (MAX(1,LWORK))
</span><span class="comment">*</span><span class="comment">          On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  LWORK   (input) INTEGER
</span><span class="comment">*</span><span class="comment">          The dimension of the array WORK. LWORK &gt;= max(1,N+M+P).
</span><span class="comment">*</span><span class="comment">          For optimum performance, LWORK &gt;= M+min(N,P)+max(N,P)*NB,
</span><span class="comment">*</span><span class="comment">          where NB is an upper bound for the optimal blocksizes for
</span><span class="comment">*</span><span class="comment">          <a name="SGEQRF.91"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</a>, <a name="SGERQF.91"></a><a href="sgerqf.f.html#SGERQF.1">SGERQF</a>, <a name="SORMQR.91"></a><a href="sormqr.f.html#SORMQR.1">SORMQR</a> and <a name="SORMRQ.91"></a><a href="sormrq.f.html#SORMRQ.1">SORMRQ</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">          If LWORK = -1, then a workspace query is assumed; the routine
</span><span class="comment">*</span><span class="comment">          only calculates the optimal size of the WORK array, returns
</span><span class="comment">*</span><span class="comment">          this value as the first entry of the WORK array, and no error
</span><span class="comment">*</span><span class="comment">          message related to LWORK is issued by <a name="XERBLA.96"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  INFO    (output) INTEGER
</span><span class="comment">*</span><span class="comment">          = 0:  successful exit.
</span><span class="comment">*</span><span class="comment">          &lt; 0:  if INFO = -i, the i-th argument had an illegal value.
</span><span class="comment">*</span><span class="comment">          = 1:  the upper triangular factor R associated with A in the
</span><span class="comment">*</span><span class="comment">                generalized QR factorization of the pair (A, B) is
</span><span class="comment">*</span><span class="comment">                singular, so that rank(A) &lt; M; the least squares
</span><span class="comment">*</span><span class="comment">                solution could not be computed.
</span><span class="comment">*</span><span class="comment">          = 2:  the bottom (N-M) by (N-M) part of the upper trapezoidal
</span><span class="comment">*</span><span class="comment">                factor T associated with B in the generalized QR
</span><span class="comment">*</span><span class="comment">                factorization of the pair (A, B) is singular, so that
</span><span class="comment">*</span><span class="comment">                rank( A B ) &lt; N; the least squares solution could not
</span><span class="comment">*</span><span class="comment">                be computed.
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">  ===================================================================
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     .. Parameters ..
</span>      REAL               ZERO, ONE
      PARAMETER          ( ZERO = 0.0E+0, ONE = 1.0E+0 )
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Local Scalars ..
</span>      LOGICAL            LQUERY
      INTEGER            I, LOPT, LWKMIN, LWKOPT, NB, NB1, NB2, NB3,
     $                   NB4, NP
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Subroutines ..
</span>      EXTERNAL           SCOPY, SGEMV, <a name="SGGQRF.123"></a><a href="sggqrf.f.html#SGGQRF.1">SGGQRF</a>, <a name="SORMQR.123"></a><a href="sormqr.f.html#SORMQR.1">SORMQR</a>, <a name="SORMRQ.123"></a><a href="sormrq.f.html#SORMRQ.1">SORMRQ</a>, <a name="STRTRS.123"></a><a href="strtrs.f.html#STRTRS.1">STRTRS</a>,
     $                   <a name="XERBLA.124"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. External Functions ..
</span>      INTEGER            <a name="ILAENV.127"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>
      EXTERNAL           <a name="ILAENV.128"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a> 
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Intrinsic Functions ..
</span>      INTRINSIC          INT, MAX, MIN
<span class="comment">*</span><span class="comment">     ..
</span><span class="comment">*</span><span class="comment">     .. Executable Statements ..
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Test the input parameters
</span><span class="comment">*</span><span class="comment">
</span>      INFO = 0
      NP = MIN( N, P )
      LQUERY = ( LWORK.EQ.-1 )
      IF( N.LT.0 ) THEN
         INFO = -1
      ELSE IF( M.LT.0 .OR. M.GT.N ) THEN
         INFO = -2
      ELSE IF( P.LT.0 .OR. P.LT.N-M ) THEN
         INFO = -3
      ELSE IF( LDA.LT.MAX( 1, N ) ) THEN
         INFO = -5
      ELSE IF( LDB.LT.MAX( 1, N ) ) THEN
         INFO = -7
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Calculate workspace
</span><span class="comment">*</span><span class="comment">
</span>      IF( INFO.EQ.0) THEN
         IF( N.EQ.0 ) THEN
            LWKMIN = 1
            LWKOPT = 1
         ELSE
            NB1 = <a name="ILAENV.159"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="SGEQRF.159"></a><a href="sgeqrf.f.html#SGEQRF.1">SGEQRF</a>'</span>, <span class="string">' '</span>, N, M, -1, -1 )
            NB2 = <a name="ILAENV.160"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="SGERQF.160"></a><a href="sgerqf.f.html#SGERQF.1">SGERQF</a>'</span>, <span class="string">' '</span>, N, M, -1, -1 )
            NB3 = <a name="ILAENV.161"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="SORMQR.161"></a><a href="sormqr.f.html#SORMQR.1">SORMQR</a>'</span>, <span class="string">' '</span>, N, M, P, -1 )
            NB4 = <a name="ILAENV.162"></a><a href="ilaenv.f.html#ILAENV.1">ILAENV</a>( 1, <span class="string">'<a name="SORMRQ.162"></a><a href="sormrq.f.html#SORMRQ.1">SORMRQ</a>'</span>, <span class="string">' '</span>, N, M, P, -1 )
            NB = MAX( NB1, NB2, NB3, NB4 )
            LWKMIN = M + N + P
            LWKOPT = M + NP + MAX( N, P )*NB
         END IF
         WORK( 1 ) = LWKOPT
<span class="comment">*</span><span class="comment">
</span>         IF( LWORK.LT.LWKMIN .AND. .NOT.LQUERY ) THEN
            INFO = -12
         END IF
      END IF
<span class="comment">*</span><span class="comment">
</span>      IF( INFO.NE.0 ) THEN
         CALL <a name="XERBLA.175"></a><a href="xerbla.f.html#XERBLA.1">XERBLA</a>( <span class="string">'<a name="SGGGLM.175"></a><a href="sggglm.f.html#SGGGLM.1">SGGGLM</a>'</span>, -INFO )
         RETURN
      ELSE IF( LQUERY ) THEN
         RETURN
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Quick return if possible
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.EQ.0 )
     $   RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Compute the GQR factorization of matrices A and B:
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">            Q'*A = ( R11 ) M,    Q'*B*Z' = ( T11   T12 ) M
</span><span class="comment">*</span><span class="comment">                   (  0  ) N-M             (  0    T22 ) N-M
</span><span class="comment">*</span><span class="comment">                      M                     M+P-N  N-M
</span><span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     where R11 and T22 are upper triangular, and Q and Z are
</span><span class="comment">*</span><span class="comment">     orthogonal.
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="SGGQRF.195"></a><a href="sggqrf.f.html#SGGQRF.1">SGGQRF</a>( N, M, P, A, LDA, WORK, B, LDB, WORK( M+1 ),
     $             WORK( M+NP+1 ), LWORK-M-NP, INFO )
      LOPT = WORK( M+NP+1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Update left-hand-side vector d = Q'*d = ( d1 ) M
</span><span class="comment">*</span><span class="comment">                                             ( d2 ) N-M
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="SORMQR.202"></a><a href="sormqr.f.html#SORMQR.1">SORMQR</a>( <span class="string">'Left'</span>, <span class="string">'Transpose'</span>, N, 1, M, A, LDA, WORK, D,
     $             MAX( 1, N ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
      LOPT = MAX( LOPT, INT( WORK( M+NP+1 ) ) )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Solve T22*y2 = d2 for y2
</span><span class="comment">*</span><span class="comment">
</span>      IF( N.GT.M ) THEN
         CALL <a name="STRTRS.209"></a><a href="strtrs.f.html#STRTRS.1">STRTRS</a>( <span class="string">'Upper'</span>, <span class="string">'No transpose'</span>, <span class="string">'Non unit'</span>, N-M, 1,
     $                B( M+1, M+P-N+1 ), LDB, D( M+1 ), N-M, INFO )
<span class="comment">*</span><span class="comment">
</span>         IF( INFO.GT.0 ) THEN
            INFO = 1
            RETURN
         END IF
<span class="comment">*</span><span class="comment">
</span>         CALL SCOPY( N-M, D( M+1 ), 1, Y( M+P-N+1 ), 1 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Set y1 = 0
</span><span class="comment">*</span><span class="comment">
</span>      DO 10 I = 1, M + P - N
         Y( I ) = ZERO
   10 CONTINUE
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Update d1 = d1 - T12*y2
</span><span class="comment">*</span><span class="comment">
</span>      CALL SGEMV( <span class="string">'No transpose'</span>, M, N-M, -ONE, B( 1, M+P-N+1 ), LDB,
     $            Y( M+P-N+1 ), 1, ONE, D, 1 )
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Solve triangular system: R11*x = d1
</span><span class="comment">*</span><span class="comment">
</span>      IF( M.GT.0 ) THEN
         CALL <a name="STRTRS.234"></a><a href="strtrs.f.html#STRTRS.1">STRTRS</a>( <span class="string">'Upper'</span>, <span class="string">'No Transpose'</span>, <span class="string">'Non unit'</span>, M, 1, A, LDA,
     $                D, M, INFO )
<span class="comment">*</span><span class="comment">
</span>         IF( INFO.GT.0 ) THEN
            INFO = 2
            RETURN
         END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">        Copy D to X
</span><span class="comment">*</span><span class="comment">
</span>         CALL SCOPY( M, D, 1, X, 1 )
      END IF
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     Backward transformation y = Z'*y
</span><span class="comment">*</span><span class="comment">
</span>      CALL <a name="SORMRQ.249"></a><a href="sormrq.f.html#SORMRQ.1">SORMRQ</a>( <span class="string">'Left'</span>, <span class="string">'Transpose'</span>, P, 1, NP,
     $             B( MAX( 1, N-P+1 ), 1 ), LDB, WORK( M+1 ), Y,
     $             MAX( 1, P ), WORK( M+NP+1 ), LWORK-M-NP, INFO )
      WORK( 1 ) = M + NP + MAX( LOPT, INT( WORK( M+NP+1 ) ) )
<span class="comment">*</span><span class="comment">
</span>      RETURN
<span class="comment">*</span><span class="comment">
</span><span class="comment">*</span><span class="comment">     End of <a name="SGGGLM.256"></a><a href="sggglm.f.html#SGGGLM.1">SGGGLM</a>
</span><span class="comment">*</span><span class="comment">
</span>      END

</pre>

 </body>
</html>
